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ABSTRACT 

The temperature of the gas in molecular clouds is a key determinant of the characteristic mass of star 
formation. Ambipolar diffusion (AD) is considered one of the most important heating mechanisms in 
weakly ionized molecular clouds. In this work, we study the AD heating rate using 2-fluid turbulence 
simulations and compare it with the overall heating rate due to turbulent dissip ation. We find that 
for observed molecular clouds, which typically have Alfven Mach numbers of ~ 1 ( Crutcher|1999 | and 
AD Reynolds numbers of ^ 20 ( McKee et al.|2010| , about 70% of the total turbulent dissipation is in 
the form of AD heating. AD has an important effect on the length scale where energy is dissipated: 
when AD heating is strong, most of the energy in the cascade is removed by ion- neutral drift, with a 
comparatively small amount of energy making it down to small scales. We derive a relation for the 
AD heating rate that describes the results of our simulations to within a factor of two. Turbulent 
dissipation, including AD heating, is generally less important that cosmic-ray heating in molecular 
clouds, although there is substantial scatter in both. 

Subject headings: ISM: kinematics and dynamics — ISM: magnetic fields — magnetic fields — magnetohydrodynamics 
(MHD) — stars:formation 



1. INTRODUCTION 

The temperature of the gas in molecular clouds is readily observable and is a key determinant of the characteristic 
mass of star formation, since the Jeans mass varies as the 3/2 power of the temperature. In molecular gas, which is 
generally shielded from far-ultraviolet radiation, the heating processes that determine the temperature are cosmic-ray 
ionization, turbulent dissipation, ambipolar diffusion, magnetic reconnect ion, and hydrody nami c compression. Cosmic 
rays are generally assumed to be the dominant heating mechanism (e.g., 'Goldsmith 2001), but Pan & Padoan (20091 
found that turbulent dissipation can give a higher heating rate in some cases. Ambipolar diffusion (AD J is the snppage 



of magnetized ions through the dominant neutral gas, and the collisions resulting from this relative motion ge nerate 
heat. AD heating is considered to be one of the most important heatins 
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significant. 

Although AD heating and turbulent dissipation have generally been considered as separate processes, in a turbulent 
medium AD heating is a component of turbulent dissipation. Turbulence is driven on large scales and cascades down 
to smaller scales. This picture is rigorously true for incompressible t urbulence and has b een shown to be consistent 
with numerical simulations for supersonic (compressible) turbulence ( Kritsuk et al.||2007 ). 

the 



AD occurs on the scale of 



the neutral-ion mean free path, which for weakly ionized gas is much greater than the neutral-neutral mean free path 
that governs viscous dissipation. As a result, energy is drained from the turbulent cascade first by AD and then by 
viscous dissipation. Magnetic reconnection also contributes to turbulent dissipation on small scales. 

In view of the importance of AD heating in turbulent media, we have re-evaluated its magnitude with higher- 
resolution simulations than were possible a decade ago. The AD heating rate per unit volume. Fad, due to this 
ion-neutral friction is 

Lad = lABPiPn (vi - v„) = jAuPiPnVa- (1) 

where pi^n and are the density and velocity of the ion and neutral components, respectively, 7ad = {o'v)/{mn + mi) 
is the ion-neutral coupling constant, and Vd is the magnitude of the ion-neutral drift velocity. The parameter (av) is 
the collision rate coefficient between ionic and neutral species, where m„ and rui are the mean neutral and ion masses, 
respectively. We are interested in the case in which the ion mass fraction, Xi = Pi/ Pin is sufficiently low that the 
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ion inertia is negligible. Solving the full two-fluid MHD equations is computationally prohibitive in this case, since 

f . —1/2 

the ion Alfven velocity va oc Xi is very large and the Courant limit on the time step correspondingly very small. 
Two approximations have been used to treat this problem: one is to use a single-fluid approximation that turns the 
induction equation into a diffusion equation; the difficulty with this approach is that the time step scales as the square 
of the grid spacing, and becomes prohibitively small at high resolution. The other appro ach is to solve the full 2-fluid 



MHD equation s with a semi-implicit approach that uses the heavy-ion approximation (P.S. Li et al. 2006 Oishi & 
Mac Low 2006), in which the ion mass fraction, x^, is raised and the coupling coefficient, 7ad, is lowered in a way 
that preserves their product, leaving the ion-neut ral drag unchanged . We hav e used this approac h to perfor m high 
resolution non-ideal MHD turbulence simulations (| P.S. Li et a r||2008| (Paper I), |McKee et al.||2010| (Paper H), |P.S Li 
(Paper IH)). It allows us to use explicit time stepping for both fluids without running afoul of the (Jl'L 



et al. 2012 (Paper IH)). 



constraint due to high ion Alfven speeds. P.S. Li et al., (2006, ) give a justification of this approach as well as a more 



detailed description of the equations and numerical tecfimques employed 

In this paper, we evaluate the importance of AD heating and compare it with the overall turbulent dissipation 
rate. In Section 2, we summarize our numerical models and the assumptions adopted. In Section 3, we report our 
simulation results and demonstrate convergence. In Section 4, we provide an analytic model that predicts the mean 
AD heating rate and show that it is consistent with the numerical results. This model gives the approximate heating 
rate in terms of observable molecular cloud properties only, so that the importance of AD drift heating can be gauged 
for real molecular clouds. In Section 5 we discuss our results a nd compare our prediction on the AD heating rate to 
the work of PZNOO, as updated in PZN12 ( Padoan et al.|[2012 ). Section 6 summarizes our conclusions. 



2. NUMERICAL MODELS AND CONVERGENCE STUDY 

The results discussed in this paper are based on the series of 512'^ AD, ideal MHD, and hydrodynamic turbulent box 
simulations investigated in Papers II and HI, plus several additional AD runs with different magnetic field strengths 
and resolutions. In all runs, the gas was isothermal, there was no gravity, and the boundaries were periodic. In the 
MHD runs, the field was initially uniform. The field strength is characterized by the plasma (3 parameter. 
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where p is the mean density in the box, Cg is the isothermal sound speed, A^a is the Alfven Mach number and M. is 
the three-dimensional rms sonic Mach number. In Papers I-III we did not distinguish the equilibrium value of /? from 
the initial value, /3o because -Brms did not differ significantly from the initial field strength, i?o, but here we do. 

To characterize the importance of AD in these simulations, we use the AD Reynolds number, which is the ratio of the 
characteristic A D timescale to the flow time, or , equivalently, the ratio of the size of the system to the characteristic 
AD length scale ( [Zweibel fc Brandenb"urg||1997[ ), 



i?AD(4) 



rms ^AD 



(3) 



Here, pi.„ is the mean density of the ions and neutrals, respectively, Iq is the size of the system (for the simulation, it 
is the size of the turbulent box), Brms is the rms magnetic fleld, and 



pv^dV 



(4) 



is the density- weighted mean-squared velocity of the system. The dynamical crossing time is t/ = ^o/'f^rms, and ^ad is 
the AD time scale. We denote the initial value of the AD Reynolds number by i?AD,o(^o)- 

The turbulence in all of the models wa s main tained at a constant rms Mach number, Al, by a flxed driving pattern 
generated using the recipe of Mac Low (1999). The turbulence was driven between wavenumbers fc = 1 ~ 2 (all 
wavenumbers are in units of 27r/£o) for a period of 3t/. We allowed the turbulence to develop for It/, and averaged all 
results over 14 data dumps spread over the subsequent 2tf unless otherwise indicated. Both the neutral and ionized 
components of the gas were driven using the same driving pattern and amplitude in order to prevent the driving from 
creating an artificial velocity difference between the two components. In the ideal run and five main AD runs, the 
turbulent box was initially threaded by a uniform magnetic field with /3o = 0.1. 

In addition to the flve main sub-Alfvenic AD models, we carried out three additional AD simulations. Two models 
had = 3 but /3o = 1 and 10 to provide information on slightly and highly super-Alfvenic turbulence. A third AD 
model driven to At = 10 at 256'^ resolution provides a high thermal Mach number model for comparison. Table 1 
summarizes the parameters of all the runs in this paper. 

The AD models are based on the assumption that the ions are a separately conserved fluid. A discussion of the 
effect of assuming ionization equilibrium instead can be found in the Appendix of Paper I and in Paper II. For our AD 
models, there is no major difference between ion conservation and ionization equilibrium for most of the turbulence 
statistics, except for the clear differences in the ion density PDF. With ionization equilibrium, the velocity power 
spectra are about 1 a steeper than for the ion conservation case. The differences in the properties of the clumps 
investigated in the above papers between the two ionization models were at the few percent level. Comparison of the 
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TABLE 1 
Model Parameters 



ModeP 


7AD 


-Rad,o(^o) 




/3o 




-Brms / Bo 


Ma 




resolution'^ 


mSph 





















0.70 


512^ 


m3c2r-l 


4 


0.12 


0.12 


0.1 


0.1 


1.00 


0.67 


0.71 


512^ 


m3c2r0 


40 


1.2 


1.2 


0.1 


0.1 


1.00 


0.67 


0.73 


512^ 


m3c2rl 


400 


12 


11.6 


0.1 


0.097 


1.02 


0.66 


0.88 


512^ 


m3c2r2 


4000 


120 


113 


0.1 


0.096 


1.03 


0.66 


0.96 


512^ 


m3c2r3 


40000 


1200 


1110 


0.1 


0.095 


1.03 


0.65 


0.92 


5123 


m3i 


oo 


oo 


oo 


0.1 


0.092 


1.04 


0.65 


0.87 


5123 


m3c2rlb0 


40 


12 


8.7 


1.0 


0.75 


1.15 


1.84 


0.86 


512^ 


m3c2rlbl 


4 


12 


6.5 


10 


5.4 


1.36 


5.07 


0.99 


512^ 


m3c2r2bl 


40 


120 


22 


10 


1.84 


2.33 


3.07 


1.01 


512^ 


ml0c2 


40 


4 


3.7 


0.1 


0.092 


1.04 


2.14 


0.77 


256^ 



Models are labeled as "mxcyrn," where x is the thermal Mach number, y = | logxioli E'nd n = log(i?AD{^)/1.2). Model "m3i" is an ideal 
MHD and "m3ph" is a pure hydrodynamic model. Model m3c2r0 is the same as model m3c2h in LMKF. Models m3c2rlbl and m3c2rlb2 
have different plasma /3o from the other models. 

Except model ml0c2, all AD models have 128^ and 256^ resolution runs for convergence study. 




log X; 

Fig. 1. — Convergence study of the volume mean AD heating rate, (Fad)) on resolution (diamonds for 128^, circles for 256^, and squares 
for 512^) and Xi- Using Xi = 0.01 is accurate enough for the heavy-ion approximation when computing (Fad)i a-nd 512^ resolution appears 
to be well-converged. The data point from the 512^ model is slightly shifted to the right and the 128^ model to the left of Xi = 0-01 to 
prevent the error bars from overlapping. 



AD heating rates with ion conservation vs. ionization equihbrium for our five main AD models, all at 256^ resolution 
shows that the differences are ^ 10 — 15%, which is about 1 a. 

We performed a convergence study of the time- and volume-averaged AD heating rate, (Fad)! in both spatial resolu- 
tion and ion-mass fraction (xi) to ensure that the results are spatially resolved and that the heavy- ion approximation 
is accurate for our chosen value of Xi = 0.01. To test the heavy- ion approximation, we used the four 256'^ models, 
mSclrO (corresponding to M — 3, Xi — 10~^ and -Rad,o(^o) = 1-2 x 10°) to m3c4r0 {xi = 10"'*), from Paper I and 
plot (Fad) vs. Xi i^i Figure [l] (Keep in mind that the physical values of Xi can be less than 10^^.) The results show 
that that using Xi = 0.01 in the heavy-ion approximation gives a value of Fad that is accurate to within 10% in this 
case. We have also run 128'^ and 256'^ simulations for the five main AD models to study the convergence behavior of 
(Fad) with spatial resolution. The values of (Fad) for the 128^ and the 512'^ simulations of model m3c2r0 are shown 
in Figure [l] for illustration. The heating rate shows the expected second-order convergence, within the uncertainties, 
and we conclude that a numerical resolution of 512^ provides a well-converged AD heating rate. We obtained similar 
results for all but the -Rad.o(^o) ~ 1200 case, which showed convergence at an order of only ~ 1.5. We will use 
Richardson extrapolation ( )Richardson||1927p to estimate the converged values of the AD heating rate in the modeling 
in Section m 



3. RESULTS FOR THE TURBULENT DISSIPATION RATE AND THE AD HEATING RATE 
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Our niain objectives are to determine the total time-averaged rate of dissipation of turbulent energy, (Tt), and that 
part of it that is due to ambipolar diffusion, (Fad), as a function of the AD Reynolds number, i?AD(^o)- We begin 
with the total time-averaged rate of dissipation of turbulent energy, which can be written as 

(r,) = 6, (5) 

where id is the typical outer length scale of the turbulence; in numerical simulations, this is the typical length scale 
on which the turbulence is driven, which in our case is £o/\/2 (fc = 1 ~ 2). In studies of incompressible turbulence, 
the dissipation rate is often normalized to the integral length scale, Ljnf We have chose n to normalize to the typical 



driving scal e, id, as has been done in some past simulations of supersonic turbulence (e.g., Mac Low|1999 Lemaster & 
Stone|2009 |) From our pure HD, ideal MHD, and all AD turbulence models at 512^, we find that (Lint)/4 = 0.85±0.15 



We have not bothered to distinguish Umis from its time average, (wrms), since the turbulent driving forces them to be 
equal to within 0.01% 

First, we discuss the turbulent dissipation rate in the HD and ideal MHD cases. Our HD model gives et 0.65, which 



is consisten t with the results of other turbulence simulations at similar resolution (e.g., Mac Low 
Stone|2009 ). We have not carried out a convergence study of our hydrodynamic result, but we note t 



19991 pemaster k 
lat tne dissipation 

rate for this simulation agrees with that for the low -Rad(^o) s imulation m3c2r-l, wh ich we verified is converged. The 



dissipation rate for incompressible HD turbulence is smaller: Kaneda et al. (2003) carried out simulations with far 
higher resolution (up to 4096^) for this case and found ~ 0.4 — 0.5. They used the integral length scale in computing 
the dissipation rate; for our pure HD model, this is almost the same as id- Our ideal MHD model gives et — 0.25, 
which is at the low end of Mac Low's (1999) results (et ~ .3 — 0.6) from his lower re solution simulations of supersonic 



MHD turbulence. In their simulations of such turbulence, Lemaster & Stone (2009) found et — 0. 4 — 0.5, about 50% 
larger than the value we find. This is probably because ot the difference m the driving method: [Lemaster fc Stone] 
( [2009j ) used a variable, sharply peaked driving pattern, while we used a fixed, flat-top driving pattern. They measured 
the dissipation time in terms of the flow time across Apk, the wavelength of the peak in the perturbation spectrum; 
this is equivalent to replacing id in Equation ^ by Apk. They found that the normalized dissipation time was very 
insensitive to changes in Apk. Studies of the energy dissipation rate in simulations of supersonic driven turbulence with 
and without magnetic fields have found that t he energy dissipation rate for ideal MHD turbulence is smaller than that 
for HD turbulence by a factor of 1.4 ~ 2 (e.g., Mac Low|1999 Lemaster fc Stone|2009 ), and our result is at the upper 



end of this range. 

The total dissipation rate is equal to the sum of the AD dissipation rate, Fad , and the other forms of dissipation, 
which we shall group together in Fother- In the simulations, Fother is a combination of numerical viscosity and numerical 
resistivity that arise from our imperfect discretization of the fluid equations. In real partially ionized, turbulent fluids 
such as molecular clouds, Fother represents energy that cascades down to very small scales and is dissipated by physical 
processes such as molecular viscosity and Ohmic dissipation. Our convergence study shows that whatever the details 
of the numerical dissipation, it occurs on sufficiently small scales in the simulation that it has a small effect on the 
mean AD heating rate. 

To facilitate comparison with the total dissipation rate, we normalize the time-averaged AD heating rate, (Fad), by 
Pi^rmj^d, so that 

(Fad)=6ad^. (6) 

We compute eAD for each of the five AD models with f3o = 0.1 and plot the results versus i?AD(^o) in Figure[2] We also 
plot the normalized total energy injection rate, et, based on the time-averaged energy driven into the box to maintain 
the gas at Mach 3. The uncertainty in e is given by the standard error of the mean, which we calculate as the standard 
deviation evaluated for a total 14 data sets in the last 2tf divided by the square root of the number of independent 
samples of e, which we estimate as 2. 

For the five AD models, et increases with decreasing i?AD(^o) from 0.25 (the same as the ideal MHD value) at the 
highest value of i?AD(^o) to 0.84 when ii!AD(^o) ~ 1; it then drops back to about 0.65 (the same as the HD value) 
as i?AD(^o) — ^ 0- The increase in et at intermediate values of i?AD(^o) is clearly due to the additional energy lost to 
AD heating. The AD heating is small in the model with the highest value of i?AD(^o) (— 1000), which is expected as 
ions and neutrals are strongly coupled. When -Rad(^o) becomes lower and the drift velocity larger, (Fad) increases, 
reaching a maximum around the models with -Rad(^o) = 11-6 and 1.2. At yet lower values of i?AD(^o), the drag 
velocity saturates at the rms velocity and the ion-neutral collision rate declines (in molecular clouds, this would most 
likely be due to a decrease in the ionization), resulting in a decreasing value of ej. Figure |3] illustrates the saturation 
of the density-weighted \vd\- As i?AD(^o) becomes small, the distribution of \vd\ almost overlaps the neutral velocity 
distribution. The volume-weighted velocity distributions also show the same trend of saturation of the magnitude of 
the drift velocity, \vd\- Figure [3[f) shows the density-weighted Vd,rms and Wn^rms of the five AD models. 

4. MODELING THE AD HEATING RATE 

Here we develop an analytic model for the AD heating rate, which depends on the ion-neutral drift velocity Vd — 
Vi — v„. The drag force on the neutrals is Fdrag — "/AuPiPn'^d, where 7ad is the ion-neutral coupling coefficient. The 
rate at which heat is generated by the drag is Fdrag • v^, so that the volume-averaged AD heating rate is (Equation 
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Fig. 2. — Normalized AD heating rates (squares) and total energy injection rates (circles) of the five AD models with /3o = 0.1 as a 
function of /?ad(^o). The energy injection rates of an ideal MHD (down triangle) model and a pure hydrodynamics (up triangle) model 
are also plotted for comparison. The -Rad(^o) of these two limiting cases are moved from oo to 10* for the ideal MHD model and from 
to 10-2 for the pure hydro model to fit on the plot. The total energy injection rates of the two limiting cases match the AD models with 
largest and smallest Rad- The drop in AD heating rate in the model with iJAD(^o) — 0.1 (m3c2r-l) is due to the saturation of the drag 
velocity in the weakly coupled ion and neutral components: the drop in the models with -Rad(^o) — 10^, 10^ is due to the reduction in 
when the ions and neutrals are well-coupled. See Section Is] for discussion. 
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FlG. 3. — Distributions of the magnitudes of the drag velocity, (solid line) and the neutral velocity, \vn\ (dashed line), of models 
(a) m3c2r-l (Rad,o(^o) = 0.12), (b) m3c2rD (ii:AD,o(^o) = 1-2), (c) m3c2rl (Had,o(^o) = 12), (d) m3c2r2 (Rad,o(^o) = 120), and (e) 
m3c2r3 (i?AD,o{^o) = 1200). (f) The normalized dispersion in the drift velocity (circles), fdrms/cs, and in the neutral velocity (squares), 
Vnrms/cs ~ as functions of -Rad(^o). 
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(1)) 



(Tad) = ( 7ADPiPnfrms ( ) ) > (7) 

~7ADPiPnW,2ms ( ) > (8) 

\ ^rms / 

where 

{vj} = {pv'a)/P (9) 

is mass-averaged, just like v'^^^ (note that (u^) is the only case in which (x) represents a mass average; otherwise it 
denotes a volume average). The second step is accurate for small and moderate values of i?AD(^d), but can be in error 
by up to a factor 3 at large i?AD(^d)- We can rewrite this in terms of r, the ratio of the flow time at the driving scale, 
— ^d/vrms, to the neutral-ion collision time, 



_ tf _ RAoi^d 

trn M 

and obtain 



4.1. ESTIMATE OF (vj) 
When the ion inertia is negligible, the equation of motion for the ions is 

lADPiPn^d = Fi, (12) 

where 

Fl--^[(VxB)xB] (13) 



47r 

is the Lorentz force. The mean-squared drag velocity is then 



1/ pWl? 



P \ {lADPtPny 



(14) 



from equation (12 1. We divide the magnetic field into a steady component and a variable one, 

B = (B)+(5B, (15) 

so that 
Now, 



AnesB ' 



(17) 



where £sb is the characteristic scale over which the fluctuating component of the field varies. If the heating were wave 
dissipation, £sb would be proportional to k~^. We now define iss by inserting this approximate expression into the 
expression for the mean-squared drag velocity. 



f>2 Sjf>2 

{vD ^ ^'"^-'•'f,^ ■ (18) 



Letting cfiB = (5i?rms/-Brms, wc find 



To evaluate we express the energy in the fluctuating field in terms of the kinetic energy as 

^^Bfms = l^PvLs, (20) 

where the factor ^ measures the deviations from equipartition (Heitsch et al. 2001; Paper III). Normalizing with 
respect to the rms field, we have 

<Pl ^ = ^Ma'. (21) 



Ambipolar Diffusion Heating in Turbulent Systems 



7 



At low Alfven Mach nu mbers, we expect that th e energy in the fluctuating field will be in equipartition with the 
associated kinetic energy ( Zweibel fc McKee|1995 ). For a uniform field, this kinetic energy is in the directions normal 



to the field, so ^ ~ |. The velocity field in the simulations is found to be approximately isotropic. As discussed above, 
this should remain valid in the presence of AD provided that the ions and neutrals are well-coupled along the field, 
which they are under the conditions we consider here. On the other hand, since i/)^ cannot exceed unity, ^ ~ Xj M.p^ 
for large values of A^a- As a result, we estimate 

2 

i- V (22) 

2 



which agrees with the results of Stone et al. (1998) to within a factor 1.2. The corresponding result for is 

(23) 



,2 ?-^a' 



The final quantity to be evaluated in Equation (19) is I^b- In the limit of large RADi^d), ^sb should have a well- 
defined value, which we term £sb, oo- We will make the assumption that isB,oo is proportional to the driving scale, 
with a constant of proportionality called a: 

i5B,ao=a£d. (24) 

The motivation for this is that in a periodic box, gradients cannot exist on scales larger than £o/2t:, and if most 
of the power is at large scales, the gradients in Isb.oo should be dominated by the largest allowed wavelengths. If 
£d ~ ^0, this implies a « l/27r. As discussed in section 4.2 below, our simulations show a ~ 0.17, which is indeed 
close to l/27r. Note that while we do find a weak dependence of £sb,oo on Ma, we are ignoring it in our analytic 
treatment, as including it does not lead to an improvement in the accuracy of our formula in the range of parameter 
space (0.7 < Ma < 5.0) that we have considered. We were unable to obtain converged results for higher equilibrium 
values of Ma with our 512^ simulations, so we were unable to confirm the 1/Ma scaling found by PZNOO 



For small values of i?AD(^o), AD will smooth the field fluctuations, so that £sb > £53,00- Hence, equation (18 1 
provides an upper limit on (v^) if £sb, 00 replaces £sb' 



2 



v-^, ^ Sa^RAMHl + IMa') ia^MA^il + lMA')r^' ^^^^ 

This equation provides a good estimate of (w^) at high values of RA'D{£d), but it is too large at moderate or low 
values of RAu{£d)- First, (w^) must be less than v'^^^; in fact, in the weak coupling limit, (w^) = | v'^^^. In Paper II, 
we introduced several different regimes of AD: Regime 1 is ideal MHD, Regime 2 is standard AD {tf > tni ^ tin), 
Regime 3 is strong AD {tni > tf > tin). Regime 4 is weak coupling {tin > tf) and Regime 5 is the hydrodynamic 
limit {tf/tin — > for fixed Ma)- In Paper II, we found that most molecular clouds are in Regime 2. Here we shall 
consider regimes 1-3, omitting regime 4, in which the heavy-ion approximation breaks down, and regime 5 in which it 
is irrelevant. In regimes 1-3, the ions are reasonably well-coupled parallel to the field, so that AD is primarily normal 
to the field. We therefore require 

{Vd)<lv^ms, (26) 

where the inequality approaches equality in the strong AD limit. We also require that the AD dissipation rate be less 
than the total dissipation rate ((Pad) < (r*)), so that 

(vj) < - vLs- (27) 

T 

In principle, we could include the factor | here also, but it does not improve the accuracy of the fit to the simulations. 

4.2. THE AD HEATING RATE 
We have now found three upper limits on (w^); one of them is a good estimate of the value of (u^) at large values 



of RAT>{£d) (Equation (25l), and one is a good estimate at small values of RAi){£d) (Equation(26)). Combining these 



two estimates and ensuring that (v^) is less than all three upper limits, we adopt 

3 T- 3 2 2 2 



""^^ - ^ ' ^ ' ^a'MA' ( 1 + Ima^ ] t\ (28) 



so that 



(rAD) ^ , ^ , — , ^ f ^ ) . (29) 
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Fig. 4. — Normalized AD heating rates of (a) five AD models (m3c2r-l to m3c2r3) with initial Mrms = 3 and (Sq = 0.1 (solid circles), (b) 
Model m3c2rlb0 with initial Alrms = 3 and /3o = 1.0 (solid down triangle), (c) Model m3c2rlbl with initial yVfrms = 3 and I3q = 10 (solid 
diamond), (d) Model ml0c2 with initial A^rms = 10 and /3o = 0.1 (solid up triangle ), a nd (e) Model m3c2r2bl with initial A^rms = 3 and 
Po = 10 (solid square). The open symbols are the predicted values from Equation ( |29| , using et = 0.55 and the rms Ma at equilibrium. 
All predictions are within a factor of two of AD heating rates measured from simulations. 



We now consider the values that the normahzed AD heating rate, ad) / {pv^msl ^d), takes in various Hmits. Varying 
one quantity at a time, we find that the normalized AD heating rate approaches for Wrms — >■ 0; for v-cms — > oo; for 
^rms — >■ 0; etT/(|et + r) for Brms oo; and for both p — >• and p — >■ oo, holding the fractional ionization constant 
in both cases. 

The careful reader will recall the comment that equation ([8) is not very accurate at large i?AD(^d)- In order to 
overcome this problem, we determine a by fitting our final result (Eq. (29 1 or pip) to our simulation results. Fitting 
the results at the largest value of i?AD(^o) that is well converged {RxD^Joj ~ 11^, we find 

a -0.17. (30) 

Our ideal MHD models at /? = 0.01, 0.1, and 1.0 show that a is almost independent of M.a- 
The AD heating rate in Equation ( 29 1 can then be rewritten as 



(r 



1 



ad; 



3Ma^ 
2RaM) 



0.043i?AD(^<i) [1 + -Ma 



(31) 



which explicitly shows that it is always less than the turbulent dissipation rate. The value of et is shown in Figure 
[2j it has a value ~ 0.65 for low values of i?AD(^o) (including the hydrodynamic case), rises to a maximum ~ 0.84 
at i?AD(^o) — 10 and then falls to about 0.3 for large values of i?AD(^o) (including the ideal MHD limit). We plot 
the normalized AD heating rates measured fro m our AD simulations versus RAT){^d) in Figure H together with the 
predicted AD heating rates from Equation ([29]) using et ~ 0.55, the mean value from all our models listed in Table 1. 
The predicted AD heating rates are all withm a factor of two of the simulation values. 
Numerically, the total turbulent dissipate rate (Equation [5| is 



3.8 X 10" 



V0.5/ 



^d, pc 



erg cm 



(32) 



where njj is the density of H nuclei, = Wrms/(10^ cm s ^), and -^d.pc = ^(i/(l pc). In estimating the heating rates in 
observed clouds, one can replace id with the observed cloud size, ^ohs- Gas in molecular clouds typically has A^a 1 
( Crutcher|p99l ) and Rad%) ~ 20 ( |McKee et al]|2010[ ). Figure |2] shows that et ~ 0.5 for this value of i?AD(M. 
For these parameters, the AD heating rate is only shghtly lower than the total dissipative heating: Equation ( [3l] ) 
implies eAD/e* — 0.8, whereas the numerical results for the model with the closest set of parameters (m3c2rl, whTch 
has i?AD,o(4) = 12) gives exo/et ~ 0.7. 



5. DISCUSSION 



How important is AD heating compared to cosmic-ray heating? Since AD heating is a substantial fraction of 
turbulent dissipation for conditions that are typical in molecular clouds, we compare with the latter. The turbulent 
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TABLE 2 

Parameters of AD simulations A1-A7 in PZN 



Name Ma P|> Rad{(-o) Fad ' 

(fiG) (km s~^) (erg cm"-^ s"-*-) (erg cm"-^ ^~^) 



Al 


83.1 


0.3 


9.9 


1.9 


4.3e5 


7.5e-30 


2.9e-33 


A2 


18.2 


1.6 


12.0 


2.3 


1.7e4 


2.7e-28 


2.7e-30 


A3 


8.2 


2.6 


9.1 


1.7 


4.6e3 


4.3e-28 


2.1e-29 


A4 


5.5 


4.5 


10.5 


2.0 


1.8e3 


1.8e-27 


1.8e-28 


A5 


4.2 


6.0 


10.8 


2.0 


1.0e3 


3.2e-27 


5.8e-28 


A6 


2.5 


7.5 


8.1 


1.5 


4.8e2 


2.7e-27 


1.3e-27 


A7 


0.7 


44.1 


12.4 


2.3 


2.4el 


2.6e-25 


1.5e-25 



"Computed for 10 K gas. 
''From Equation (31). 



line width is obse rved to increase with scale as cr = v^^as/V'^ = o'pc(^pc/2)'''^, where typically dpc — 0.7 km s ^ (McKee 
fc Ostriker||2007 1 . The turbulent heating rate in molecular clouds is then 

r,.2.4xlO-(iL)(;j^^)"r.„f3t »gc„-s-. (33) 

By comparison, the cosmic-ray ionization rate in molecular clouds with colum n densities A^h S 2 x 10^ ^ cm~^, which 
has recently been revised upwards, is about 3.5 x 10~^^ s~^ (H2 molecule)"^ ( Indriolo fc McCaIIl|2012 ). Glassgold et 



al. (20121 find that the heating per ionization in molecular clouds is 13 eV, so this ionization rate corresponds to a 
heating rate of Fcr — 3.6 x lO^^^'nn erg cm""^ s~^. Under normal circumstances then, turbulent dissipation, including 
AD heating, is not competitive with cosmic-ray heating in molecular clouds. There are significant variations in both 
the linewidth-size relation and the cosmic-ray heating rate, however, so turbulent dissipation can dominate in some 
regions. For very large molecular clouds {£d ^ 100 pc), the two heating rates become comparable, but it should be 
borne in mind that the heating due to turbulent dissipation is spatially localized: A single shock wave at a velocity 



Wrms that extends across the area of the cloud. A, could account for the bulk of the heating due to turbulent dissipation 
in the cloud, since the volume-averaged heating rate is ^pyf^^^A/V = ^pv'^^^/i ~ Ft, where V is the volume of the 
cloud and i its size. 

Finally, we compare our results on the AD heating rate in turbulent systems to those of PZNOO, who examined 
this question previously. Note that PZNOO originally contained numerical errors that caused them to overestimate the 
mean AD heating rate. Therefore, we compare against their corrected numbers in PZN12 (when it is not important 
to distinguish between the two papers, we shall simply refer to PZN). We focus on their runs A1-A7, which they used 
to determine an expression for the mean AD heating rate. For convenience, we have summarized the key physical 
parameters of these runs in Table 2. 

PZN characterized the strength of AD in their simulations through the parameter a, which is related to our i?AD(^o) 
by 

Ml N 

^Mio) = ^-, (34) 

where iV is their numerical resolution. All the runs in Table 2 have N = 128 and a — 0.21, but the range in i?AD(^o) is 
very large. Note, however, that PZN have no runs with i?AD(^o) ^ 24, so we cannot compare to them in that regime. 
We can, however, check for consistency at higher i?AD(^o)- 
From their simulations, PZN12 infer the mean AD heating rate in a turbulent molecular cloud for the case in which 

Observe that we have expressed the density in terms of the number of hydrogen nuclei, rather than the number of 
neutral particles, as they did. PZNOO estimate that Fad oc a"-^ over the range 0.11 < a < 0.74, but this is questionable 
since a depends on the numerical resolution, whereas the AD heating rate does not. 
To compare their results with ours, we compute the mean heating rate twice for eac h run listed in Table 2, once 



using their result (Equation |35[ ), which we label FpzN, and once using our Equation (31 1, which we simply label Fad- 
The results are listed in Table 2. We find that the agreement is fairly good (within a factor ~ 5) for A^a ^ 4 and 
-Rad(-^o) < 1000, which is the regime we explored directly. However, for their runs at high Ma and high RadC^o): the 
disagreement is much worse, and by run Al, we are lower by about 3 orders of magnitude. 

What accounts for the difference? While PZN find that the magnetic length scale is proportional to 1/A^a, we find 
it to be only weakly dependent on A^a in the range we were able to consider (A^a < 5.0), and treat it as constant in 
our analytic theory. Although we could improve our agreement with PZN by inserting a 1/A^a dependence (actually, 
we would need a 1/(1 + Ma) dependence to be consistent with our results at low A^a) into the third term in the 



denominator in Equation (31 1, we have chosen not to do so, since our simulation results do not provide evidence for 
this dependence, at least over the limited range of A^a we could explore. Note that run Al has i?AD(^o) = 4.3 x 10^, 
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meaning that this run was very close to the ideal MHD limit. The AD heating rate in this regime is only a negligible 
fraction of the total turbulent heating. The regime in which AD heating is significant corresponds to the last few 
runs in Table 2, which is where our agreement is fairly good, especially considering the differences in our numerical 
approaches — e.g., we use a two-fluid treatment to their one, we assume ion conservation and they that the ionization 
is determined by cosmic rays, and our turbulence driving used a fixed velocity pattern whereas theirs was random in 
time and was applied as a source term in the moment um equation. Finall y, we note that real molecular clouds are 
observed to lie mainly in the moderate i?AD(^o) regime McKee et al. (2010). Thus, the agreement is reasonably good 
in the regimes of the most physical and astrophysical interest. 

6. CONCLUSIONS 

We have discussed a number of the effects of ambipolar diffusion (AD) on weakly magnetized molecular clouds in 
Papers II and III, based on high resolution, two-fluid MHD turbulence simulations using the heavy-ion approximation. 
The strength of AD can be measured with the AD Reynolds number, i?AD(^o) (Equation [s]). In the limit of low 
-Rad(^o)j one recovers the hydrodynamic limit, while in the limit of high i?AD(^o) one recovers ideal MHD. Molecular 
clouds are observed to have intermediate values of -Rad(^o)i ran ging from 3 to about 70 in Crutcher's (1999) sample of 
molecular clouds with measured field strengths ( McKee et al.|2010 ). In this paper, we focus on the heating due to the 
friction between ions and neutrals inside weakly ionized, turbulent molecular clouds over a wide range of conditions. 
We compute the AD heating rates directly from our two-fluid model using Equation ([l]). Our conclusions are: 

1. We find that the AD heating rate is a significant fraction of the overall turbulent heating rate in the range of 
-Rad(^o) = 1 ^ 100, provided Ma is not large. As noted above, this range of i?AD (^o) encompasse s observed 



(| Crutcher|p99| . Our AD 
ion is m the form of AD 



molecular clouds; furthermore, observed molecular clouds typically have ^ 1 
turbulence simulations show that in this regime, up to 70-80% of the total dissipaf 
heating. AD heating gives a moderate increase in the total turbulent dissipation rate in molecular clouds with 
typical values of i?AD(^o)- At smaller values (i?AD(^o) ^ 1), the AD heating rate falls off rapidly as a result of 
infrequent collisions between ions and neutrals. 

2. Heating due to turbulent dissipation, including AD heating, is generally less than cosmic-ray heating in molecular 
clouds, although there is substantial scatter in both rates. 

3. AD significantly affects the length scale at which turbulent energy is dissipated. When AD is weak, either because 
the two fluids are too weakly coupled to cause significant momentum exchange (i?AD(^o) ^ 1) or because they 
are so well-coupled that no drift develops (-Rad(^o) ^ 100) (FigurelSl), almost all of the energy can cascade down 
to the viscous and/or resistive scales. On the other hand, when AT) heating is strong (i?AD(-^o) ~ 1 — 100, the 
regime most significant for molecular clouds), we find that most of the energy in the cascade can be removed by 
ion- neutral drift, with < | of the turbulent energy making it down to small scales. The scale at which turbulent 



energy is converted to heat can affect the temperature distribution of interstellar gas (Pan & Padoan 2009j); 
however, we defer discussion of this effect to future work. 



4. We derive a relation (Equation 31 ) for the AD heating rate based on the global physical properties of the system. 
The predicted ratio of the AlJneating rate to the total is accurate to within a factor of two for all our AD 
models. The relation is useful in predicting the AD heating rate inside molecular clouds from their observed 
properties. 
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